library(BiocManager)
source("rnaseq_utils.R")
library(DESeq2)
library(BiocParallel)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(tidyr)
library(edgeR)
library(sys)
# remove the outlier
rnaCounts <- rnaCounts %>% select(-`4GEXDARK4`)
sampleAnnotation <- sampleAnnotation[colnames(rnaCounts),] 
riboCounts <- riboCounts %>% select(-`4GEXDARK4`)
sampleAnnotation2 <- sampleAnnotation2[colnames(riboCounts),] 
rnaCounts = data.frame(rnaCounts)
rownames(sampleAnnotation) = colnames(rnaCounts)
sampleAnnotation$group <- ifelse(substr(sampleAnnotation$group, 1, 1) %in% c("1", "4"),
                                 paste0("X", sampleAnnotation$group), 
                                 sampleAnnotation$group) 

RNA-Seq

EdgeR

group = factor(sampleAnnotation$group)
y = DGEList(counts = rnaCounts, 
        group = group,
        genes = rownames(rnaCounts))
y$samples
##                  group lib.size norm.factors
## X14BENDDAY2 X14BENDDAY  9022678            1
## X14BENDDAY4 X14BENDDAY  8646052            1
## X14BEXDARK2 X14BEXDARK  6505229            1
## X14BEXDARK3 X14BEXDARK  8245123            1
## X14BEXDARK4 X14BEXDARK  7509975            1
## X4GENDDAY2   X4GENDDAY  8258548            1
## X4GENDDAY3   X4GENDDAY  8113002            1
## X4GENDDAY4   X4GENDDAY  7815864            1
## X4GEXDARK2   X4GEXDARK  6772977            1
## X4GEXDARK3   X4GEXDARK  7203593            1
## COLENDDAY3   COLENDDAY  7986754            1
## COLENDDAY5   COLENDDAY  9449171            1
## COLEXDARK2   COLEXDARK  6979346            1
## COLEXDARK3   COLEXDARK  7910292            1
## COLEXDARK4   COLEXDARK  7288743            1
design <- model.matrix(~0+group)
colnames(design) <- levels(group)
keep <- filterByExpr(y, design)
y <- y[keep, , keep.lib.sizes=FALSE]
AveLogCPM <- aveLogCPM(y)
hist(AveLogCPM)

y <- calcNormFactors(y)
y$samples
##                  group lib.size norm.factors
## X14BENDDAY2 X14BENDDAY  9020237    0.9304865
## X14BENDDAY4 X14BENDDAY  8643156    0.9110121
## X14BEXDARK2 X14BEXDARK  6503769    0.8708110
## X14BEXDARK3 X14BEXDARK  8242249    0.9027839
## X14BEXDARK4 X14BEXDARK  7508091    0.8972249
## X4GENDDAY2   X4GENDDAY  8255380    1.1210043
## X4GENDDAY3   X4GENDDAY  8109195    1.1550759
## X4GENDDAY4   X4GENDDAY  7811866    1.1826089
## X4GEXDARK2   X4GEXDARK  6770891    0.9697061
## X4GEXDARK3   X4GEXDARK  7201373    0.9719922
## COLENDDAY3   COLENDDAY  7982688    1.1459087
## COLENDDAY5   COLENDDAY  9445634    1.0863755
## COLEXDARK2   COLEXDARK  6977274    0.9617534
## COLEXDARK3   COLEXDARK  7907452    0.9781382
## COLEXDARK4   COLEXDARK  7286107    0.9894661
for (i in 1:ncol(y)){
  plotMD(y, column=i)
  abline(h=0, col="red", lty=2, lwd=2)
}

contrast between day and dark; wild-type

start_time <- Sys.time()

group = factor(sampleAnnotation$group)
y = DGEList(counts = rnaCounts, 
        group = group,
        genes = rownames(rnaCounts))
#y$samples

design <- model.matrix(~0+group)
colnames(design) <- levels(group)
keep <- filterByExpr(y, design)
y <- y[keep, , keep.lib.sizes=FALSE]
AveLogCPM <- aveLogCPM(y)
#hist(AveLogCPM)
y <- calcNormFactors(y)
#y$samples

y <- estimateDisp(y, design, robust=TRUE) 
#y$samples
#plotBCV(y) 

fit <- glmQLFit(y, design, robust=TRUE) 
#plotQLDisp(fit)
#summary(fit$df.prior)

day_vs_dark <- makeContrasts(COLENDDAY - COLEXDARK, levels = design) 
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient:  1*COLENDDAY -1*COLEXDARK 
##               genes     logFC unshrunk.logFC   logCPM       PValue          FDR
## AT5G20630 AT5G20630  5.535501       5.535991 9.724466 7.830552e-15 1.302143e-10
## AT1G69570 AT1G69570 -4.169657      -4.176277 4.964722 5.986341e-14 4.120473e-10
## AT3G09600 AT3G09600 -6.186618      -6.203965 5.354701 8.471142e-14 4.120473e-10
## AT5G05580 AT5G05580  3.479189       3.482670 5.151173 9.911536e-14 4.120473e-10
## AT1G44446 AT1G44446 -3.324531      -3.325194 7.499148 1.560527e-13 4.843196e-10
## AT5G45820 AT5G45820 -6.180542      -6.196619 5.706870 1.747500e-13 4.843196e-10
## AT1G07180 AT1G07180 -4.592792      -4.600768 5.580593 2.696248e-13 6.405131e-10
## AT3G27690 AT3G27690 -4.871827      -4.872516 8.909493 3.335744e-13 6.933760e-10
## AT5G64840 AT5G64840 -3.751972      -3.752280 8.916823 9.196219e-13 1.649972e-09
## AT3G59060 AT3G59060 -3.863711      -3.864380 7.955721 1.012013e-12 1.649972e-09
is.de <- decideTestsDGE(tr)
summary(is.de)
##        1*COLENDDAY -1*COLEXDARK
## Down                       1959
## NotSig                    12389
## Up                         2281
plotMD(tr, status=is.de)

output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 4998
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])

end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 7.519514 secs
start_time <- Sys.time()

DESeqDataSet = DESeqDataSetFromMatrix(
    countData = rnaCounts,
    colData = sampleAnnotation,
    design = ~ group
)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
DESeqDataSet = DESeq(
  DESeqDataSet, 
  parallel=FALSE
  ) 
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLENDDAY", "COLEXDARK"))
# indexes that all have a value 
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 9513
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 951.3
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 4998
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 18.03931 secs

contrast between day and dark; 14G mutant

start_time <- Sys.time()

day_vs_dark <- makeContrasts(X14BENDDAY - X14BEXDARK, levels = design) 

tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient:  1*X14BENDDAY -1*X14BEXDARK 
##               genes     logFC unshrunk.logFC   logCPM       PValue          FDR
## AT5G45820 AT5G45820 -9.429530      -9.603817 5.706870 5.558380e-14 6.183772e-10
## AT3G27690 AT3G27690 -5.451716      -5.452763 8.909493 8.843228e-14 6.183772e-10
## AT3G47500 AT3G47500 -6.387256      -6.397688 6.374169 1.156934e-13 6.183772e-10
## AT2G05070 AT2G05070 -6.702382      -6.703587 9.838101 1.534871e-13 6.183772e-10
## AT1G69570 AT1G69570 -3.999885      -4.006312 4.964722 1.859334e-13 6.183772e-10
## AT5G57660 AT5G57660 -5.172456      -5.173075 9.479774 3.123604e-13 8.657069e-10
## AT1G18330 AT1G18330 -4.796010      -4.803761 5.709382 4.024328e-13 9.560079e-10
## AT1G44446 AT1G44446 -3.090320      -3.090875 7.499148 5.304359e-13 1.043978e-09
## AT3G59060 AT3G59060 -4.096335      -4.097397 7.955721 5.650253e-13 1.043978e-09
## AT1G07180 AT1G07180 -3.544550      -3.546441 5.580593 2.491965e-12 4.143888e-09
is.de <- decideTestsDGE(tr)
summary(is.de)
##        1*X14BENDDAY -1*X14BEXDARK
## Down                         1880
## NotSig                      12443
## Up                           2306
plotMD(tr, status=is.de)

output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 4954
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])

end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.256883 secs
start_time <- Sys.time()

DESeq_Results <- results(DESeqDataSet, contrast = c("group", "X14BENDDAY", "X14BEXDARK"))
# indexes that all have a value 
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 9667
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 966.7
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 4954
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.470524 secs

contrast between day and dark; 4G

start_time <- Sys.time()

day_vs_dark <- makeContrasts(X4GENDDAY - X4GEXDARK, levels = design) 

tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient:  1*X4GENDDAY -1*X4GEXDARK 
##               genes     logFC unshrunk.logFC   logCPM       PValue          FDR
## AT1G44446 AT1G44446 -3.953663      -3.954660 7.499148 1.525969e-15 2.537534e-11
## AT5G45820 AT5G45820 -6.767170      -6.793169 5.706870 1.180308e-14 5.976616e-11
## AT1G69570 AT1G69570 -4.069463      -4.075807 4.964722 1.453358e-14 5.976616e-11
## AT3G09600 AT3G09600 -5.679975      -5.691259 5.354701 1.502934e-14 5.976616e-11
## AT1G18330 AT1G18330 -4.701840      -4.706428 5.709382 1.797046e-14 5.976616e-11
## AT3G27690 AT3G27690 -5.016305      -5.017302 8.909493 2.333974e-14 6.468610e-11
## AT3G59060 AT3G59060 -4.007284      -4.007909 7.955721 7.243031e-14 1.720634e-10
## AT3G47500 AT3G47500 -5.279691      -5.286332 6.374169 8.482539e-14 1.763202e-10
## AT5G57660 AT5G57660 -4.701959      -4.702394 9.479774 1.015748e-13 1.876763e-10
## AT2G05070 AT2G05070 -5.289302      -5.289925 9.838101 2.110553e-13 3.266227e-10
is.de <- decideTestsDGE(tr)
summary(is.de)
##        1*X4GENDDAY -1*X4GEXDARK
## Down                       1720
## NotSig                    12872
## Up                         2037
plotMD(tr, status=is.de)

output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 4476
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])

end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.352772 secs
start_time <- Sys.time()

DESeq_Results <- results(DESeqDataSet, contrast = c("group", "X4GENDDAY", "X4GEXDARK"))
# indexes that all have a value 
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 9236
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 923.6
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 4476
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.394359 secs

contrast between genotypes, day: COL and 14B

start_time <- Sys.time()

day_vs_dark <- makeContrasts(COLENDDAY - X14BENDDAY, levels = design) 

tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient:  1*COLENDDAY -1*X14BENDDAY 
##               genes     logFC unshrunk.logFC   logCPM       PValue          FDR
## AT5G57870 AT5G57870  2.749030       2.749739 7.240689 3.247423e-10 4.619118e-06
## AT5G13210 AT5G13210 -3.824220      -3.838854 3.198047 5.555497e-10 4.619118e-06
## AT5G44610 AT5G44610  2.040782       2.042746 5.585858 1.396575e-08 7.741215e-05
## AT2G18193 AT2G18193 -3.343663      -3.348822 3.958703 5.256919e-08 2.185433e-04
## AT5G60660 AT5G60660  2.159697       2.162756 4.918551 7.241155e-08 2.408263e-04
## AT4G34131 AT4G34131 -1.938129      -1.938531 6.332505 1.258730e-07 3.488570e-04
## AT3G01345 AT3G01345 -3.602524      -3.631246 5.944690 1.924598e-07 4.572020e-04
## AT4G17340 AT4G17340  1.757094       1.757392 7.348337 2.432282e-07 5.055802e-04
## AT4G25210 AT4G25210 -1.541807      -1.542088 6.480979 2.947914e-07 5.221489e-04
## AT3G03920 AT3G03920 -1.887654      -1.887918 6.613104 3.139990e-07 5.221489e-04
is.de <- decideTestsDGE(tr)
summary(is.de)
##        1*COLENDDAY -1*X14BENDDAY
## Down                         284
## NotSig                     16218
## Up                           127
plotMD(tr, status=is.de)

output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 623
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])

end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.322701 secs
start_time <- Sys.time()

DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLENDDAY", "X14BENDDAY"))
# indexes that all have a value 
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 4011
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 401.1
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 623
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.027226 secs

contrast between genotypes, night: COL and 14B

start_time <- Sys.time()

day_vs_dark <- makeContrasts(COLEXDARK - X14BEXDARK, levels = design) 

tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient:  1*COLEXDARK -1*X14BEXDARK 
##               genes     logFC unshrunk.logFC   logCPM       PValue          FDR
## AT5G13210 AT5G13210 -3.699243      -3.711515 3.198047 6.465266e-11 1.075109e-06
## AT5G57870 AT5G57870  2.260281       2.260651 7.240689 4.047196e-10 3.365041e-06
## AT5G60660 AT5G60660  2.054796       2.055994 4.918551 2.406564e-09 1.333958e-05
## AT5G14250 AT5G14250 -1.632070      -1.632924 4.981151 3.681860e-09 1.530641e-05
## AT5G66170 AT5G66170  2.316835       2.319725 3.877768 4.723212e-09 1.570846e-05
## AT4G12545 AT4G12545  4.856443       4.889670 3.900695 6.141792e-09 1.702198e-05
## AT4G22490 AT4G22490  2.334567       2.336895 5.016727 1.271862e-08 2.880894e-05
## AT3G21720 AT3G21720 -3.113993      -3.117116 4.479165 1.385961e-08 2.880894e-05
## AT3G11710 AT3G11710 -1.509525      -1.509708 7.236135 2.197153e-08 3.541515e-05
## AT5G43290 AT5G43290 -3.499606      -3.512768 2.504409 2.236003e-08 3.541515e-05
is.de <- decideTestsDGE(tr)
summary(is.de)
##        1*COLEXDARK -1*X14BEXDARK
## Down                         289
## NotSig                     15992
## Up                           348
plotMD(tr, status=is.de)

output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 929
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])

end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.830286 secs
start_time <- Sys.time()

DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLEXDARK", "X14BEXDARK"))
# indexes that all have a value 
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 5578
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 557.8
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 929
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.36721 secs

contrast between genotypes, day: COL and 4G

start_time <- Sys.time()

day_vs_dark <- makeContrasts(COLENDDAY - X4GENDDAY, levels = design) 

tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient:  1*COLENDDAY -1*X4GENDDAY 
##               genes     logFC unshrunk.logFC   logCPM       PValue          FDR
## AT3G01345 AT3G01345 -8.047863      -8.079019 5.944690 3.056858e-14 5.083250e-10
## AT4G15110 AT4G15110  2.878915       2.880247 5.730227 1.877357e-12 1.560928e-08
## AT3G60240 AT3G60240  1.843232       1.843563 7.002984 7.683085e-09 4.258734e-05
## AT1G71030 AT1G71030  2.093515       2.094985 7.033848 2.007862e-07 8.347183e-04
## AT1G80840 AT1G80840  2.046024       2.047474 4.407891 7.730338e-07 2.570956e-03
## AT4G23810 AT4G23810  2.186970       2.190009 3.329671 2.564529e-06 7.107593e-03
## AT1G68840 AT1G68840  1.340636       1.340815 7.553711 6.326135e-06 1.493487e-02
## AT5G03230 AT5G03230  1.488172       1.489980 6.029975 7.184977e-06 1.493487e-02
## AT4G38400 AT4G38400 -1.721493      -1.725036 3.439118 9.668246e-06 1.786370e-02
## AT5G61600 AT5G61600  1.581549       1.582239 5.208417 1.652684e-05 2.748248e-02
is.de <- decideTestsDGE(tr)
summary(is.de)
##        1*COLENDDAY -1*X4GENDDAY
## Down                          3
## NotSig                    16615
## Up                           11
plotMD(tr, status=is.de)

output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 15
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])

end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.248233 secs
start_time <- Sys.time()

DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLENDDAY", "X4GENDDAY"))
# indexes that all have a value 
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 352
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 35.2
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 15
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.037442 secs

contrast between genotypes, night: COL and 4G

start_time <- Sys.time()

day_vs_dark <- makeContrasts(COLEXDARK - X4GEXDARK, levels = design) 

tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient:  1*COLEXDARK -1*X4GEXDARK 
##               genes     logFC unshrunk.logFC   logCPM       PValue          FDR
## AT3G01345 AT3G01345 -8.971074      -9.050976 5.944690 6.666672e-15 1.108601e-10
## AT3G60240 AT3G60240  1.632912       1.633235 7.002984 1.923606e-07 1.599382e-03
## AT3G10340 AT3G10340  1.638733       1.639960 4.514862 3.316951e-05 1.838586e-01
## AT5G23010 AT5G23010  1.423243       1.423636 7.184509 1.770548e-04 6.469499e-01
## AT3G21720 AT3G21720 -1.718978      -1.721435 4.479165 1.945246e-04 6.469499e-01
## AT5G13210 AT5G13210 -1.626017      -1.635016 3.198047 3.212706e-04 8.904015e-01
## AT2G40330 AT2G40330 -1.978164      -1.983003 3.000185 4.745570e-04 9.779622e-01
## AT2G19970 AT2G19970 -1.405319      -1.406831 3.186981 5.036444e-04 9.779622e-01
## AT5G24240 AT5G24240 -2.511893      -2.532129 1.051515 5.902314e-04 9.779622e-01
## AT2G45360 AT2G45360 -1.608946      -1.615313 1.778893 6.644364e-04 9.779622e-01
is.de <- decideTestsDGE(tr)
summary(is.de)
##        1*COLEXDARK -1*X4GEXDARK
## Down                          1
## NotSig                    16627
## Up                            1
plotMD(tr, status=is.de)

output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 2
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])

end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.288385 secs
start_time <- Sys.time()

DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLEXDARK", "X4GEXDARK"))
# indexes that all have a value 
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 355
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 35.5
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 2
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.407555 secs

contrast between genotypes, day: 14B and 4G

start_time <- Sys.time()

day_vs_dark <- makeContrasts(X14BENDDAY - X4GENDDAY, levels = design) 

tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient:  1*X14BENDDAY -1*X4GENDDAY 
##               genes     logFC unshrunk.logFC   logCPM       PValue          FDR
## AT4G15110 AT4G15110  2.774434       2.775750 5.730227 4.013744e-12 3.956349e-08
## AT3G01345 AT3G01345 -4.445340      -4.447772 5.944690 4.758373e-12 3.956349e-08
## AT5G57870 AT5G57870 -2.774680      -2.775392 7.240689 1.331251e-10 7.379125e-07
## AT3G60240 AT3G60240  1.830524       1.830854 7.002984 8.904515e-09 3.701829e-05
## AT5G13210 AT5G13210  2.497301       2.502468 3.198047 2.013228e-08 5.781913e-05
## AT3G03920 AT3G03920  2.005648       2.005943 6.613104 2.217572e-08 5.781913e-05
## AT4G12545 AT4G12545 -4.246320      -4.261516 3.900695 2.433904e-08 5.781913e-05
## AT4G34131 AT4G34131  1.860616       1.860990 6.332505 4.188871e-08 8.707092e-05
## AT3G07230 AT3G07230  2.189971       2.190546 5.937643 6.721483e-08 1.098499e-04
## AT2G18193 AT2G18193  2.796600       2.799935 3.958703 7.443747e-08 1.098499e-04
is.de <- decideTestsDGE(tr)
summary(is.de)
##        1*X14BENDDAY -1*X4GENDDAY
## Down                         221
## NotSig                     15910
## Up                           498
plotMD(tr, status=is.de)

output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 1013
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])

end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.812748 secs
start_time <- Sys.time()

DESeq_Results <- results(DESeqDataSet, contrast = c("group", "X14BENDDAY", "X4GENDDAY"))
# indexes that all have a value 
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 5629
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 562.9
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 1013
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.369587 secs

contrast between genotypes, night: 14B and 4G

start_time <- Sys.time()

day_vs_dark <- makeContrasts(X14BEXDARK - X4GEXDARK, levels = design) 

tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient:  1*X14BEXDARK -1*X4GEXDARK 
##               genes     logFC unshrunk.logFC   logCPM       PValue          FDR
## AT3G01345 AT3G01345 -4.947873      -4.952555 5.944690 3.957324e-13 6.580635e-09
## AT4G22490 AT4G22490 -2.904688      -2.907204 5.016727 5.414356e-10 4.501766e-06
## AT5G57870 AT5G57870 -2.164870      -2.165234 7.240689 2.806398e-09 1.555587e-05
## AT5G13180 AT5G13180 -2.014076      -2.014778 5.870159 1.080500e-08 4.491907e-05
## AT5G14250 AT5G14250  1.697728       1.698642 4.981151 1.387473e-08 4.614457e-05
## AT3G60240 AT3G60240  1.862560       1.862905 7.002984 2.101704e-08 5.437957e-05
## AT4G12545 AT4G12545 -4.713775      -4.746878 3.900695 2.289115e-08 5.437957e-05
## AT2G31280 AT2G31280 -1.620274      -1.621100 4.937200 4.808899e-08 9.912306e-05
## AT2G45670 AT2G45670 -1.618117      -1.618455 6.162460 5.364770e-08 9.912306e-05
## AT3G07230 AT3G07230  2.405400       2.406257 5.937643 6.682140e-08 1.111173e-04
is.de <- decideTestsDGE(tr)
summary(is.de)
##        1*X14BEXDARK -1*X4GEXDARK
## Down                         289
## NotSig                     16062
## Up                           278
plotMD(tr, status=is.de)

output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 861
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])

end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.285785 secs
start_time <- Sys.time()

DESeq_Results <- results(DESeqDataSet, contrast = c("group", "X14BEXDARK", "X4GEXDARK"))
# indexes that all have a value 
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 5456
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 545.6
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 861
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.365066 secs

Ribo-Seq

riboCounts = data.frame(riboCounts)
rownames(sampleAnnotation2) = colnames(riboCounts)
sampleAnnotation2$group <- ifelse(substr(sampleAnnotation2$group, 1, 1) %in% c("1", "4"),
                                 paste0("X", sampleAnnotation2$group), 
                                 sampleAnnotation2$group) 
group = factor(sampleAnnotation2$group)
y = DGEList(counts = riboCounts, 
        group = group,
        genes = rownames(riboCounts))
y$samples
##                  group lib.size norm.factors
## X14BENDDAY2 X14BENDDAY   994611            1
## X14BENDDAY4 X14BENDDAY  1128982            1
## X14BEXDARK2 X14BEXDARK   407864            1
## X14BEXDARK3 X14BEXDARK   564737            1
## X14BEXDARK4 X14BEXDARK   332326            1
## X4GENDDAY2   X4GENDDAY   890314            1
## X4GENDDAY3   X4GENDDAY   925586            1
## X4GENDDAY4   X4GENDDAY  1258788            1
## X4GEXDARK2   X4GEXDARK  1151971            1
## X4GEXDARK3   X4GEXDARK   521701            1
## COLENDDAY3   COLENDDAY   234275            1
## COLENDDAY5   COLENDDAY   123746            1
## COLEXDARK2   COLEXDARK   556287            1
## COLEXDARK3   COLEXDARK  2133574            1
## COLEXDARK4   COLEXDARK   350392            1
design <- model.matrix(~0+group)
colnames(design) <- levels(group)
keep <- filterByExpr(y, design)
y <- y[keep, , keep.lib.sizes=FALSE]
AveLogCPM <- aveLogCPM(y)
hist(AveLogCPM)

y <- calcNormFactors(y)
y$samples
##                  group lib.size norm.factors
## X14BENDDAY2 X14BENDDAY   966881    1.0673652
## X14BENDDAY4 X14BENDDAY  1097093    1.0193377
## X14BEXDARK2 X14BEXDARK   398728    0.9595619
## X14BEXDARK3 X14BEXDARK   551606    0.9632180
## X14BEXDARK4 X14BEXDARK   325221    0.9245671
## X4GENDDAY2   X4GENDDAY   861062    1.1586712
## X4GENDDAY3   X4GENDDAY   900301    0.9996082
## X4GENDDAY4   X4GENDDAY  1220563    1.0979268
## X4GEXDARK2   X4GEXDARK  1121272    0.9923085
## X4GEXDARK3   X4GEXDARK   509437    0.8782312
## COLENDDAY3   COLENDDAY   228596    0.9600032
## COLENDDAY5   COLENDDAY   120166    1.1027311
## COLEXDARK2   COLEXDARK   542330    0.9738623
## COLEXDARK3   COLEXDARK  2073598    1.0485569
## COLEXDARK4   COLEXDARK   341882    0.8978050
for (i in 1:ncol(y)){
  plotMD(y, column=i)
  abline(h=0, col="red", lty=2, lwd=2)
}

contrast between day and dark; wild-type

start_time <- Sys.time()

group = factor(sampleAnnotation2$group)
y = DGEList(counts = riboCounts, 
        group = group,
        genes = rownames(riboCounts))
#y$samples

design <- model.matrix(~0+group)
colnames(design) <- levels(group)
keep <- filterByExpr(y, design)
y <- y[keep, , keep.lib.sizes=FALSE]
AveLogCPM <- aveLogCPM(y)
#hist(AveLogCPM)
y <- calcNormFactors(y)
#y$samples

y <- estimateDisp(y, design, robust=TRUE) 
#y$samples
#plotBCV(y) 

fit <- glmQLFit(y, design, robust=TRUE) 
#plotQLDisp(fit)
#summary(fit$df.prior)

day_vs_dark <- makeContrasts(COLENDDAY - COLEXDARK, levels = design) 
tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient:  1*COLENDDAY -1*COLEXDARK 
##               genes      logFC unshrunk.logFC    logCPM       PValue
## AT2G21660 AT2G21660   5.208113   5.216850e+00  8.733563 2.094121e-13
## AT3G63160 AT3G63160   4.949438   4.965324e+00  6.930860 2.415943e-13
## AT3G27690 AT3G27690  -5.388810  -5.396124e+00  9.575708 1.573315e-12
## AT5G20630 AT5G20630   4.556578   4.557859e+00 10.254881 3.134203e-12
## AT3G46970 AT3G46970   5.008629   5.031453e+00  7.569230 3.773531e-12
## AT4G27310 AT4G27310 -11.000459  -1.442695e+08  7.738613 1.020064e-11
## AT5G23240 AT5G23240   6.777450   6.993705e+00  5.901661 1.901113e-11
## AT2G28000 AT2G28000   3.552737   3.555919e+00  8.735217 6.338472e-11
## AT3G45300 AT3G45300  -4.940535  -4.945190e+00  9.635867 6.619262e-11
## AT4G36850 AT4G36850  -5.037495  -5.042810e+00  9.567315 6.718235e-11
##                    FDR
## AT2G21660 1.336016e-09
## AT3G63160 1.336016e-09
## AT3G27690 5.800287e-09
## AT5G20630 8.347051e-09
## AT3G46970 8.347051e-09
## AT4G27310 1.880318e-08
## AT5G23240 3.003758e-08
## AT2G28000 7.430367e-08
## AT3G45300 7.430367e-08
## AT4G36850 7.430367e-08
is.de <- decideTestsDGE(tr)
summary(is.de)
##        1*COLENDDAY -1*COLEXDARK
## Down                        903
## NotSig                     9306
## Up                          851
plotMD(tr, status=is.de)

output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 2356
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])

end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 4.851323 secs
start_time <- Sys.time()

DESeqDataSet = DESeqDataSetFromMatrix(
    countData = riboCounts,
    colData = sampleAnnotation2,
    design = ~ group
)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
DESeqDataSet = DESeq(
  DESeqDataSet, 
  parallel=FALSE
  ) 
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLENDDAY", "COLEXDARK"))
# indexes that all have a value 
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 2977
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 297.7
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 2356
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 17.4248 secs

contrast between day and dark; 14B mutant

start_time <- Sys.time()

day_vs_dark <- makeContrasts(X14BENDDAY - X14BEXDARK, levels = design) 

tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient:  1*X14BENDDAY -1*X14BEXDARK 
##               genes      logFC unshrunk.logFC   logCPM       PValue
## AT4G27310 AT4G27310  -7.611322  -7.731697e+00 7.738613 3.531205e-15
## AT3G47500 AT3G47500  -6.982655  -7.079459e+00 7.313586 2.064323e-14
## AT3G12320 AT3G12320  -7.755058  -7.913509e+00 7.591271 3.447189e-14
## AT3G27690 AT3G27690  -5.296898  -5.302185e+00 9.575708 4.741226e-14
## AT5G06980 AT5G06980 -10.868484  -1.442695e+08 7.519672 4.432385e-13
## AT2G46830 AT2G46830  -8.865830  -9.297847e+00 7.467552 8.195362e-13
## AT3G46970 AT3G46970   5.257768   5.276175e+00 7.569230 1.012280e-12
## AT2G05070 AT2G05070  -5.979803  -5.988081e+00 9.535313 1.865140e-12
## AT5G57660 AT5G57660  -4.473787  -4.477452e+00 9.302224 8.375795e-12
## AT1G06570 AT1G06570  -3.705266  -3.708559e+00 8.484350 4.465806e-11
##                    FDR
## AT4G27310 3.905513e-11
## AT3G47500 1.141571e-10
## AT3G12320 1.270864e-10
## AT3G27690 1.310949e-10
## AT5G06980 9.804435e-10
## AT2G46830 1.510678e-09
## AT3G46970 1.599402e-09
## AT2G05070 2.578556e-09
## AT5G57660 1.029292e-08
## AT1G06570 4.939182e-08
is.de <- decideTestsDGE(tr)
summary(is.de)
##        1*X14BENDDAY -1*X14BEXDARK
## Down                          870
## NotSig                       9258
## Up                            932
plotMD(tr, status=is.de)

output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 2414
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])

end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.117298 secs
start_time <- Sys.time()

DESeq_Results <- results(DESeqDataSet, contrast = c("group", "X14BENDDAY", "X14BEXDARK"))
# indexes that all have a value 
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 3385
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 338.5
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 2414
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.478439 secs

contrast between day and dark; 4G

start_time <- Sys.time()

day_vs_dark <- makeContrasts(X4GENDDAY - X4GEXDARK, levels = design) 

tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient:  1*X4GENDDAY -1*X4GEXDARK 
##               genes      logFC unshrunk.logFC   logCPM       PValue
## AT4G27310 AT4G27310  -7.701849  -7.806090e+00 7.738613 7.837052e-17
## AT3G12320 AT3G12320  -9.723506  -1.033340e+01 7.591271 3.010892e-16
## AT3G47500 AT3G47500  -7.416789  -7.560539e+00 7.313586 8.625411e-16
## AT5G06980 AT5G06980 -11.222042  -1.442695e+08 7.519672 1.246396e-14
## AT2G46830 AT2G46830  -8.710055  -9.046606e+00 7.467552 2.857466e-14
## AT3G27690 AT3G27690  -4.579005  -4.584166e+00 9.575708 4.131001e-14
## AT1G19530 AT1G19530  -7.344015  -7.425529e+00 7.135410 7.960830e-14
## AT3G16150 AT3G16150  -5.325967  -5.360825e+00 6.352284 4.201473e-13
## AT2G31380 AT2G31380  -5.048744  -5.071614e+00 7.224844 4.328159e-13
## AT1G07040 AT1G07040  -3.911955  -3.920789e+00 7.418560 4.564919e-13
##                    FDR
## AT4G27310 8.667779e-13
## AT3G12320 1.665024e-12
## AT3G47500 3.179902e-12
## AT5G06980 3.446284e-11
## AT2G46830 6.320715e-11
## AT3G27690 7.614812e-11
## AT1G19530 1.257811e-10
## AT3G16150 5.048800e-10
## AT2G31380 5.048800e-10
## AT1G07040 5.048800e-10
is.de <- decideTestsDGE(tr)
summary(is.de)
##        1*X4GENDDAY -1*X4GEXDARK
## Down                       1235
## NotSig                     8588
## Up                         1237
plotMD(tr, status=is.de)

output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 3210
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])

end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.635423 secs
start_time <- Sys.time()

DESeq_Results <- results(DESeqDataSet, contrast = c("group", "X4GENDDAY", "X4GEXDARK"))
# indexes that all have a value 
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 4242
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 424.2
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 3210
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.375332 secs

contrast between genotypes, day: COL and 14B

start_time <- Sys.time()

day_vs_dark <- makeContrasts(COLENDDAY - X14BENDDAY, levels = design) 

tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient:  1*COLENDDAY -1*X14BENDDAY 
##               genes     logFC unshrunk.logFC    logCPM       PValue        FDR
## AT5G52310 AT5G52310  2.298437       2.302264  6.242966 3.167691e-06 0.03503466
## AT2G24050 AT2G24050  3.958974       4.001360  5.655656 6.942280e-06 0.03839081
## AT5G11330 AT5G11330  3.469210       3.506837  5.021713 1.431804e-05 0.05278583
## AT5G57870 AT5G57870  2.035819       2.039778  6.949040 3.201044e-05 0.08365051
## AT1G29395 AT1G29395  3.998408       4.026022  5.913189 3.781668e-05 0.08365051
## AT4G36390 AT4G36390 -2.341185      -2.346701  6.177771 7.153424e-05 0.13186144
## AT2G30570 AT2G30570  1.862147       1.862420 10.538731 9.603520e-05 0.13473429
## AT1G51400 AT1G51400  1.689335       1.690210  8.926355 1.112608e-04 0.13473429
## AT1G54410 AT1G54410  1.722246       1.722848  9.235126 1.189421e-04 0.13473429
## AT2G15970 AT2G15970  1.837798       1.838640  8.355813 1.218212e-04 0.13473429
is.de <- decideTestsDGE(tr)
summary(is.de)
##        1*COLENDDAY -1*X14BENDDAY
## Down                           0
## NotSig                     11058
## Up                             2
plotMD(tr, status=is.de)

output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 5
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])

end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.042515 secs
start_time <- Sys.time()

DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLENDDAY", "X14BENDDAY"))
# indexes that all have a value 
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 295
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 29.5
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 5
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.609078 secs

contrast between genotypes, night: COL and 14B

start_time <- Sys.time()

day_vs_dark <- makeContrasts(COLEXDARK - X14BEXDARK, levels = design) 

tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient:  1*COLEXDARK -1*X14BEXDARK 
##               genes     logFC unshrunk.logFC   logCPM       PValue          FDR
## AT4G14130 AT4G14130  3.469643       3.472552 7.859450 1.945077e-09 2.151255e-05
## AT1G55490 AT1G55490 -2.728129      -2.731292 8.973621 3.370916e-08 1.251508e-04
## AT3G56240 AT3G56240  2.332565       2.334391 8.437367 3.924470e-08 1.251508e-04
## AT5G39610 AT5G39610  2.744453       2.751394 6.325387 4.526248e-08 1.251508e-04
## AT1G12900 AT1G12900 -2.433954      -2.434493 9.985552 1.216752e-07 2.691454e-04
## AT3G53460 AT3G53460 -2.755251      -2.759374 8.854840 2.348539e-07 4.329140e-04
## AT3G21720 AT3G21720 -4.232050      -4.291644 4.526621 2.880863e-07 4.551763e-04
## AT4G15390 AT4G15390  3.144799       3.165702 5.863056 5.068291e-07 7.006913e-04
## AT1G74880 AT1G74880 -3.140100      -3.176221 5.634535 9.189849e-07 9.218939e-04
## AT3G46010 AT3G46010  2.266082       2.273641 6.647010 9.443118e-07 9.218939e-04
is.de <- decideTestsDGE(tr)
summary(is.de)
##        1*COLEXDARK -1*X14BEXDARK
## Down                         193
## NotSig                     10710
## Up                           157
plotMD(tr, status=is.de)

output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 589
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])

end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.000142 secs
start_time <- Sys.time()

DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLEXDARK", "X14BEXDARK"))
# indexes that all have a value 
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 1786
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 178.6
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 589
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.342794 secs

contrast between genotypes, day: COL and 4G

start_time <- Sys.time()

day_vs_dark <- makeContrasts(COLENDDAY - X4GENDDAY, levels = design) 

tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient:  1*COLENDDAY -1*X4GENDDAY 
##               genes     logFC unshrunk.logFC   logCPM       PValue         FDR
## AT3G60240 AT3G60240  2.189418   2.192475e+00 7.754499 3.920502e-07 0.004336075
## AT4G36850 AT4G36850 -2.817004  -2.821714e+00 9.567315 3.956657e-06 0.021880315
## AT3G15356 AT3G15356 -2.593780  -2.600656e+00 7.013375 6.057575e-05 0.223322582
## AT5G20150 AT5G20150  2.810670   2.831369e+00 4.183695 1.308052e-04 0.301615767
## AT1G02640 AT1G02640 -1.951093  -1.954489e+00 8.492358 1.363543e-04 0.301615767
## AT2G17230 AT2G17230 -1.966528  -1.970363e+00 6.563663 2.130801e-04 0.344802627
## AT2G15080 AT2G15080 -7.394377  -1.442695e+08 4.578464 2.182295e-04 0.344802627
## AT3G45970 AT3G45970 -2.298422  -2.304392e+00 6.588248 2.583441e-04 0.357160767
## AT1G22770 AT1G22770  2.720781   2.751204e+00 4.236109 3.218990e-04 0.395578152
## AT1G15430 AT1G15430 -7.127060  -1.442695e+08 4.353627 5.250017e-04 0.548527832
is.de <- decideTestsDGE(tr)
summary(is.de)
##        1*COLENDDAY -1*X4GENDDAY
## Down                          1
## NotSig                    11058
## Up                            1
plotMD(tr, status=is.de)

output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 2
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])

end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.044354 secs
start_time <- Sys.time()

DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLENDDAY", "X4GENDDAY"))
# indexes that all have a value 
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 91
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 9.1
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 2
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.016833 secs

contrast between genotypes, night: COL and 4G

start_time <- Sys.time()

day_vs_dark <- makeContrasts(COLEXDARK - X4GEXDARK, levels = design) 

tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient:  1*COLEXDARK -1*X4GEXDARK 
##               genes     logFC unshrunk.logFC   logCPM       PValue        FDR
## AT3G60240 AT3G60240  1.862816       1.864929 7.754499 5.438739e-06 0.06015245
## AT4G20390 AT4G20390  3.323862       3.403996 4.126823 1.316746e-04 0.60463238
## AT3G27160 AT3G27160  1.680439       1.684490 7.616775 2.323116e-04 0.60463238
## AT2G40330 AT2G40330 -2.192904      -2.206413 5.029629 2.573583e-04 0.60463238
## AT3G25730 AT3G25730 -3.094740      -3.195454 3.589990 2.733419e-04 0.60463238
## AT2G43920 AT2G43920 -1.822146      -1.831592 4.344406 4.399068e-04 0.75017639
## AT5G09220 AT5G09220 -2.151239      -2.169488 4.732271 4.747952e-04 0.75017639
## AT5G23010 AT5G23010  1.666882       1.671034 7.404208 1.001170e-03 1.00000000
## AT2G23170 AT2G23170 -1.651932      -1.655450 5.954169 1.103876e-03 1.00000000
## AT4G14130 AT4G14130  1.487309       1.487828 7.859450 1.195258e-03 1.00000000
is.de <- decideTestsDGE(tr)
summary(is.de)
##        1*COLEXDARK -1*X4GEXDARK
## Down                          0
## NotSig                    11060
## Up                            0
plotMD(tr, status=is.de)

output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 1
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])

end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.103256 secs
start_time <- Sys.time()

DESeq_Results <- results(DESeqDataSet, contrast = c("group", "COLEXDARK", "X4GEXDARK"))
# indexes that all have a value 
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 25
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 2.5
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 1
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.344416 secs

contrast between genotypes, day: 14B and 4G

start_time <- Sys.time()

day_vs_dark <- makeContrasts(X14BENDDAY - X4GENDDAY, levels = design) 

tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient:  1*X14BENDDAY -1*X4GENDDAY 
##               genes     logFC unshrunk.logFC   logCPM       PValue         FDR
## AT3G60240 AT3G60240  2.121398       2.124404 7.754499 1.508249e-07 0.001409122
## AT5G57870 AT5G57870 -2.368533      -2.372729 6.949040 3.711322e-07 0.001409122
## AT2G24050 AT2G24050 -3.950367      -3.992707 5.655656 4.412224e-07 0.001409122
## AT5G17710 AT5G17710  2.269856       2.273933 6.150627 5.096283e-07 0.001409122
## AT2G18193 AT2G18193  3.057497       3.101369 4.288943 1.230352e-06 0.002721539
## AT5G36910 AT5G36910 -2.650604      -2.662379 7.378530 2.130043e-06 0.003926379
## AT4G15110 AT4G15110  5.222394       5.324235 4.877607 4.138304e-06 0.006515614
## AT5G11330 AT5G11330 -3.048942      -3.085444 5.021713 4.712922e-06 0.006515614
## AT2G22360 AT2G22360  1.838789       1.845057 5.169158 2.214719e-05 0.027216437
## AT5G08610 AT5G08610  2.370609       2.381078 4.951388 2.869953e-05 0.030705817
is.de <- decideTestsDGE(tr)
summary(is.de)
##        1*X14BENDDAY -1*X4GENDDAY
## Down                           7
## NotSig                     11043
## Up                            10
plotMD(tr, status=is.de)

output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 40
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])

end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.072463 secs
start_time <- Sys.time()

DESeq_Results <- results(DESeqDataSet, contrast = c("group", "X14BENDDAY", "X4GENDDAY"))
# indexes that all have a value 
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 502
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 50.2
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 40
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.390091 secs

contrast between genotypes, night: 14B and 4G

start_time <- Sys.time()

day_vs_dark <- makeContrasts(X14BEXDARK - X4GEXDARK, levels = design) 

tr <- glmTreat(fit, contrast=day_vs_dark, lfc=log2(1.5))
topTags(tr)
## Coefficient:  1*X14BEXDARK -1*X4GEXDARK 
##               genes     logFC unshrunk.logFC    logCPM       PValue
## AT3G56240 AT3G56240 -2.770708      -2.772654  8.437367 3.163044e-09
## AT3G12780 AT3G12780  2.651389       2.651975 10.702835 3.073265e-08
## AT3G27160 AT3G27160  2.761626       2.766593  7.616775 1.052553e-07
## AT1G55490 AT1G55490  2.874402       2.877830  8.973621 1.329417e-07
## AT2G43030 AT2G43030  2.204186       2.205813  8.520207 4.697147e-07
## AT5G39610 AT5G39610 -2.537727      -2.544460  6.325387 5.685509e-07
## AT3G51920 AT3G51920 -1.991011      -1.992451  7.715316 5.728854e-07
## AT3G46010 AT3G46010 -2.438612      -2.446364  6.647010 6.317070e-07
## AT3G24430 AT3G24430  3.012059       3.028821  6.325834 7.374642e-07
## AT2G41680 AT2G41680  2.096484       2.099069  7.449179 1.107221e-06
##                    FDR
## AT3G56240 3.498327e-05
## AT3G12780 1.699516e-04
## AT3G27160 3.675838e-04
## AT1G55490 3.675838e-04
## AT2G43030 8.733350e-04
## AT5G39610 8.733350e-04
## AT3G51920 8.733350e-04
## AT3G46010 8.733350e-04
## AT3G24430 9.062616e-04
## AT2G41680 1.054867e-03
is.de <- decideTestsDGE(tr)
summary(is.de)
##        1*X14BEXDARK -1*X4GEXDARK
## Down                         133
## NotSig                     10790
## Up                           137
plotMD(tr, status=is.de)

output = as.data.frame(topTags(tr, n = 10000))
sum(output$FDR < 0.1)
## [1] 459
sig_genes_edge <- rownames(output[which(output$FDR < 0.1),])

end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.07013 secs
start_time <- Sys.time()

DESeq_Results <- results(DESeqDataSet, contrast = c("group", "X14BEXDARK", "X4GEXDARK"))
# indexes that all have a value 
clean_DESeq_padj <- which(!is.na(DESeq_Results$padj))
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1)
## [1] 1425
sum(DESeq_Results[clean_DESeq_padj, "padj"] <= 0.1) * 0.1
## [1] 142.5
sig_genes_DESeq <- rownames(DESeq_Results[which(clean_DESeq_padj <= 0.1)])
length(intersect(sig_genes_DESeq, sig_genes_edge))
## [1] 459
end_time <- Sys.time()
print(end_time - start_time)
## Time difference of 1.362078 secs